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SUMMARY 


A generalized subsonic unsteady aerodynamic kernel function, valid for both 
growing and decaying oscillatory motions, is developed and applied in a modified 
flutter analysis computer program to solve for boundaries of constant damping ratio 
as well as the flutter boundary. Results are given for the variation of generalized 
aerodynamic forces with the damping ratio. Some comparisons are made with an alter- 
native method of obtaining generalized forces based on rational function approxima- 
tions of simple harmonic forces. For decaying motion, instances were observed of 
unexpected looping and spiraling of generalized aerodynamic forces calculated in the 
complex plane at high reduced frequencies. Similar spiraling behavior was found for 
two-dimensional flow for which convergence of the results with respect to downwash 
collocation was proven. Rates of change of damping ratios with respect to dynamic 
pressure near flutter are substantially lower from the generalized-kernel-function 
calculations than from the conventional velocity-damping (V-g) calculation. For the 
DAST ARW-1 1 , calculated values of both damping ratios and frequencies agreed with the 
in-flight experimental values for Mach numbers approaching the flutter condition. 

The aerodynamic forces from the rational function approximation used in control 
theory for s-plane analysis agreed fairly well with kernel-function results except 
for strongly damped motion at combinations of high (subsonic) Mach number and reduced 
frequency. 


INTRODUCTION 

Considerable flutter analysis has been accomplished using the subsonic kernel 
function that originated from the lifting-surface-theory work of H. G. Kussner (1940) 
that is based on the linear theory of potential flow. Watkins, Runyan, and Wools ton 

(1955) cast the kernel function in a form amenable to automated computation. The 

function was developed by assuming simple harmonic (i.e., constant-amplitude) motion 

that had continued for an infinite time. Thus, in a flutter analysis the resulting 

aerodynamic forces are valid only at the flutter boundary, but are not strictly valid 
for growing or decaying motion at speeds above and below a flutter boundary. 

Numerous analyses have attempted to forecast the near-flutter and approach-to- 
flutter behavior of lifting surfaces. Such efforts began at least as early as the 
velocity-damping (V-g) solutions of Smilg and Wasserman (1942) and continue today in 
the rational function approximation ( RFA) described by Sevart (1975) and applied, for 
example, by Abel (1979). 

The present report generalizes the subsonic kernel function into the s-plane for 
arbitrary values of the complex reduced frequency. This kernel function is then 
a PPlied to analyze the near-flutter behavior of several configurations using the 
equations of dynamic equilibrium of Cunningham (1978). The calculations were carried 
out by a modified version of the flutter analysis computer programs described by 
Desmarais and Bennett (1978). An appendix describes the evaluation of the kernel 
function by an accurate series approximation for which the computation is economical. 


wing 


1 Drones for Aerodynamic and Structural Testing (DAST), 
1 (ARW-1). 


aeroelastic research 
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The present results are compared with both the near-flutter results from the V-g 
method and the results from RFA aerodynamics in control theory* 


SYMBOLS 
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AC 


p j 


speed of sound 

generalized aerodynamic force, equation (5) 

coefficients of the rational function approximation (RFA), equation (6) 
root semichord, often reference length 
li f t-cur ve s lope 

2 

lifting pressure coefficient for unit amplitude of mode j, Ap./(pV /2) 


f 

g 

hi(x,y) 

i 

lm( ) 
k 

k c 

K(M,k,x Q , 


cyclic frequency of motion 

mechanical hysteresis structural damping coefficient of mode i 

root-finding increment added to 

mode shape of mode i 

imaginary unit, \f-T 

imaginary part of ( ) 

Wb o 

reduced frequency, and 

complex reduced frequency, k(1 + i£) 

y Q ) generalized kernel function for real k and for k replaced by 
k(1 + iO 


K(M,k,x ,y ) kernel function with singular factor removed, equation (3) 
o o 

i reference length, optionally equal to b Q 

m^ generalized mass for mode i 

M Mach number of undisturbed stream 

N number of modes in analysis 

Ap(x,y,t) lifting pressure distribution 
A Pj = 3(Ap)/3( qj /Jl) 
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q^(t) generalized coordinate of mode i 

rational function approximation for any A^_., equation (6) 
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Re( ) 

real part of ( ) 

RFA 

rational function approximation 
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= ia)( 1 + ic), 
motion e s ; 


coefficient of t for general exponential time-varying 
also the Laplace transform variable 
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= sZ/V 

planform area of lifting surface 
time 

= C-x Q + m r)/(3 2 |v o I) 

speed of undisturbed flow 


w(x,y,t) downwash distribution/ positive with z 


x,y,z 


orthogonal, right-handed coordinates nondimensionalized by 
downstream, y positive on the right-hand ha If -span, z 


Z, x positive 
positive up 


a 

3 

3 
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= x - 5 
= y - n 

angle of attack, pitch displacement 

= \fl - M 2 

aerodynamic lag parameters in equation (6) 

logarithmic decrement of oscillatory motion, positive for decay 
damping ratio of control theory, cos 0 
damping ratio used herein, defined in equation (9) 
dummy variable for y 

angle from negative real s-axis to s of motion (see sketch on p, 5) 

mass ratio, wing mass divided by the mass of air in a truncated cylindrical 
cone that just encloses the wing planform 

dummy variable for x 
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p 


air density 


a) circular frequency of motion, radians/sec 

0 ), natural frequency of mode i 

w base or reference frequency, equations (4) 

o 

u) natural frequency of torsion mode 

ot 


ANALYSIS 

Downwash Integral Equation and Kernel Function 

The phenomenon studied here is that of a thin lifting surface oscillating normal 
to the direction of a compressible flow. The linear, small disturbance equations of 
potential flow along with the associated boundary conditions are used to analyze this 
phenomenon • 

Based on the pioneer work of Kussner (1940), the integral equation that relates 
downwash and lifting pressure difference on a planar lifting surface can be put in 
the form 


-w ( x , y , t ) 
V 


— / / Ap(5,n,t) K(M,k, x ,y ) d(£5) dUn) 
4TtpV 2 S ° ° 


(1 ) 


where K(M,k,x 0 ,y 0 ) is the kernel function relating the downwash produced at x,y 
due to unit lifting pressure at £,n. 


Watkins, Runyan, and Woolston (1955) expressed the subsonic planar (z = 0) 
kernel function as 


K(M,k,x ,y ) = lim -<exp(-ikx 

o o ^ „ „ Z \ 4 

z+0 8 (&z) 


o' J 


° exp[ik(x - M\[x 2 + r 2 )/g 2 ] 


\ .2 ; 

U + r 


dX> (2) 


where X, the streamwise variable of integration for x - is positive downstream. 
The integral in equation (2) sums up the effects on the downwash at x,y at the 
present instant of all pressure disturbances which originated at over all past 

time • 


Carrying out the differentiation of equation (2), 
grating by parts changes the form of the planar kernel 


taking the limit, 
function to 


and inte- 
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where 


K(M,k,x ,y ) = exp(-ikx ) 
o o o 




exp(-iku 1 |y o | ) 



a 



exp(-ik|y o |u) 



( 3 ) 


The calculation of K is described in the appendix# 

Consider now the generalization from simple harmonic (constant-amplitude) motion 
with time variation exp iu)t (real a)) to growing and decaying oscillatory motion 
with time variation exp st (complex s). Such motions have come to be termed 
"s-plane motions," where in the complex s-plane, s = iw is the positive imaginary 
axis. As indicated in the accompanying sketch, for growing motion the real part of 
s is positive and for decaying motion the real part of s is negative. 



The generalization to growing and decaying motion from harmonic motion is accom- 
plished by replacing the frequency to by oj (1 + i£) , where £ is a motion damping 
ratio that is positive for decay. The parameter 5 = (5/2ir, where 6 is the 
logarithmic decrement of motion, as described, for example, in Cunningham (1978). 

The sketch shows that this definition of the damping ratio is equivalent to 
g = -Re(s)/Im(s) = cot 0, in contrast to the damping ratio used in control theory, 

£ = -Re ( s ) / | s | = cos 0 = i/'T + g 2 . The selection of the cotangent function to 
define the damping ratio simplifies the modifications required in the flutter analy- 
sis computer program described below. For small damping ratios the two definitions 
are asymptotic. For instance, with 0 = 70°, they differ by only 6 percent. 

The integral in equation (3) exists for growing and for constant-amplitude 
motion, but for decaying motion the integral is improper. This problem is surmounted 
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as follows* For both growing and constant-amplitude motion, the integrals defining 
the kernel function are computable. These integrals are recognized as representa- 
tions of analytic functions, namely, Bessel and Struve functions. These analytic 
functions can then be used to evaluate the quantity represented by the integrals for 
decaying motion even though the integrals are nonconvergent . This use of analytic 
continuation is described, for example, by Carrier, Krook, and Pearson (1966) and 
cited by Edwards (1977). The evaluation of the improper integral for decaying motion 
is described in the appendix. 

For the present report the FAST (Flutter Analysis System) computer program of 
Desmarais and Bennett (1978) was modified to incorporate the generalized kernel func- 
tion. The resulting aerodynamic forces are those corresponding to the user-selected 
The 12-term exponential series approximation D12.1 (see appendix) of Desmarais 
(1982) was used, in general. Several checks to verify accuracy were made with the 
24-term series D24.2. As shown in the appendix, approximation D12.1 is accurate for 
growing and decaying motion in the s-plane extending at least 45° on both sides of 
the imaginary s-axis that represents harmonic motion. 


Equations of Dynamic Equilibrium 

In equation (22) of Cunningham (1978), the equations of equilibrium are given 
for growing and decaying motion. An equilibrium condition exists when the complex 
structural forces are balanced by the complex unsteady aerodynamic forces. Solutions 
can be obtained by an application of the familiar V-g root-finding technique. In 
this method aerodynamic forces are calculated for an assumed complex reduced fre- 
quency, and an eigenvalue problem is solved to determine the roots of the equilibrium 
equations. In general, none of the roots match the assumed reduced frequency and a 
matched condition must be determined by iteration. This procedure was implemented by 
modifying the solution procedure used for the traditional V-g root-finding technique 
in the FAST computer program of Desmarais and Bennett (1978). 


Equation (5) of Desmarais and Bennett (1978) expresses the equilibrium condition 
at a flutter boundary, with motion neither growing nor decaying, for which solutions 
are obtained in the FAST program. This equation and program can be generalized to 
solve equation (22) of Cunningham (1978) for boundaries of selected £ * 0, (as well 
as for £ = 0). This generalization is done by replacing oi and k by w(1 + i£) 
and k ( 1 + i£) throughout the analysis. The resulting equations can be put in the 
form 
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0 ) \2 , 2 ^2 
o ) k (1 + ig) 
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where the generalized aerodynamic force element 
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is a function of M, k, 5 , and the planform. The eigenvalues (one for each elastic 
mode ) 


ft 



+ ig) 


are sought at values of k for which the structural damping g - 0. Traditionally, 
solutions with g $ 0 , which are obtained in the iterative process of finding a 
matched flutter point, have been used as an indication of the approach to flutter. 

For the present study the only changes in the equilibrium equations (4) from those of 
Desmarais and Bennett (1978) are the factor of (1 + i?) in the second term and the 

A. . , which are now functions of c. 

13 


A Rational Function Approximation (RFA) of Unsteady Aerodynamic Forces 

Sevart (1975), Roger (1977), and Abel (1979) described a mathematical technique 
used in control theory for approximating the unsteady aerodynamic forces in the 
general s-plane based on the forces for simple harmonic motion. The technique is 
outlined as follows: 

1 • Beginning with simple harmonic motion, each individual generalized-force 
matrix element a. of equations (4) is calculated for a series of values of ik 

(note the distinction between the modal index i (in A. . ) and the imaginary unit 
i = (in ik) ) . 13 

2. Using temporarily here the notation of Abel (1979), namely, Q(ik) to repre- 
sent any individual A. .(ik), the variation of the aerodynamic forces is approximated 

by a rational function (RFA) with real coefficients of the form 


a y (ik)A 

Q(ik) - A q + ikA 1 + (ik) Z A 2 + 2^ Ik ' 4 . -A— < 6 > 

m=3 m-2 


A 

3. Using known values of Q at a sequence of at least four values of k, a 
least-square-error solution for the real coefficients Aq to A 6 is calculated. 
(The resulting approximations can be compared with the exact (input) values of k 
and also calculated for intermediate values of k.) 


4. Substitute the complex quantity s£/V for ik in equation ( 6 ) to obtain 


A s 
m 


Q(sVV) - A Q + (s£/V)A 1 + (s£/V) A,, + / , ^ 

m=5 m -2 


(7) 


where 


s = (tf(~£ + i) = ica( 1 + i^) 


( 8 ) 
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and the damping ratio (positive for decaying motion) is 


5 = lr = cot e (9) 


In the present report all the results from the RFA were calculated using 3 m _2 as 
0.2 , 0 »4 ^ 0.6 , and 0 • 8 . 


RESULTS AND DISCUSSION 

Three flutter analyses were performed using the FAST program modified for grow- 
ing and decaying motion as described above. The planforms analyzed, depicted in 
figure 1 , are as follows: (1) a clipped-tip delta wing, (2) an aspect-ratio-five 

rectangular-planform model, and (3) a transport-type supercritical wing. Results 
of the present analysis, using the generalized kernel function, are compared with 
results of two approximate methods for calculating motion-damping ratios. One of 
these is the traditional V-g calculation of Smilg and Wasserman (1942). The other 
is the RFA aerodynamic-force method described in the preceding section. In all 
cases b Q was used for the reference length £. 


Clipped-Tip Delta-Wing Flutter Model 

The clipped-tip delta-wing flutter model analyzed is that of the sample case 
of Desmarais and Bennett (1978), and the planform and aeroelastic parameters are 
given therein. The same model was analyzed at supersonic Mach numbers in 
Cunningham (1978). Figure 2 is a plot of speed index versus density for M = 0.8. 
Kernel-function results are shown at constant damping ratios £ for harmonic 
motion (? =0), for growing motion (£ = -0.02), and for decaying motion 
(£ = 0.02 and 0.05), all with modal damping coefficients g^ = 0. For comparison 
the results of the V-g method are given by the dashed curves. In the V-g method the 

aerodynamic forces for £ = 0 are used, and the value of g^ (the same for all 

modes i) corresponds to damping ratio -2£. The dashed curves were therefore com- 
puted with the coefficients g^_ = -0.04 and -0.10 to forecast decaying motion and 
with g^ = 0.04 to forecast growing motion. 

The figure includes a possible wind-tunnel operating curve for fixed Mach num- 
ber M and speed of sound a g • The intersections of the operating curve with the 
various boundaries are projected down to an inset plot of £ and -g^/2 versus 

density. To avoid clutter, only five of eight projected points are indicated by the 

short-dashed lines. This plot reveals for this case that along the illustrated tun- 
nel operating curve near the flutter condition, the gradient of damping ratio with 
respect to increasing density is only about half as steep for the correct aerodynam- 
ics (£ variable) as it is for simple harmonic aerodynamics (£ = 0); thus the correct 
aerodynamics predict a less abrupt approach to flutter. 

Figure 3 shows curves of the speed index for constant damping ratio versus Mach 
number for M ranging from 0.5 to 0.9. The values of speed index from figure 2 
for p = 0.002378 slug/ft 3 are included. At M = 0.5, there is very little differ- 
ence between the solid-curve and the dashed-curve predictions of speed index for the 
two damping ratios analyzed, which are above and below the flutter boundary. As M 
increases toward 0.9, the two predictions become more different, with the present 
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variable-^ aerodynamics predicting a lower rate of change of damping ratio versus 
speed index at a given Mach number than does the V-g analysis with £ = 0 and 
variable g^ • 


Aspect-Ratio-Five Rectangular-Planform Model 

Aerodynamic-f orce results .- Doggett, Rainey, and Morgan (1959) tested this model 
at several Mach numbers. Yates et al. (1982) reported recent flutter analyses via 
two subsonic aerodynamic programs and compared the analytical and experimental 
results in figures 14 to 16 therein. One of the programs used was the unmodified 
FAST program of Desmarais and Bennett (1978). Three bending and two torsion modes of 
a uniform cantilever beam were utilized. 

With the present modified FAST program, aerodynamic forces were calculated for 
three decay ratios, namely, C = 0 and ±0.1, and for a range of reduced frequencies 
k from 0 to 0.228. (The flutter experiments had values of k ranging from 0.047 
to 0.102.) For one of the experimental Mach numbers, M = 0.904, two of the general- 
ized forces A and A are plotted in figure 4. The kernel-function results are 
12 22 

indicated by "K" in the figure. Figure 4(a) shows A ^ t descriptively termed the 

"weighted twisting moment due to first- tors ion-mode vibration." The curve for har- 
monic motion, £ = 0, reflects the effects induced by the past history of constant- 
amplitude motion, and the moment lags in phase behind the torsional deflection. The 
negative imaginary part indicates that this moment of itself acts to damp out the 
first torsion component of motion. The positive real part indicates that the moment 
^22 acts to i ncrease * or diverge, the torsional displacement, and thus to decrease 

the first torsional frequency as the dynamic pressure increases. The curve for 
£ = 0.1 reflects the effects induced by a past history of exponentially decaying 
motion, while the curve for £ = -0.1 reflects the effects induced by a past his- 
tory of growing motion. The three vectors from the origin to the points for 
k = 0.114 show that the moment a^ 2 lags the torsional displacement a little more 

for decaying motion and a little less for growing motion than it does for constant- 
amplitude motion. Thus to the extent that the wing torsional motion component influ- 
ences the overall flutter motion, the variation of the lag of the moment A with 

22 

£ contributes in the same sense as that is, decaying motion is accompanied by an 

increased decay-producing moment, and growing motion by a decreased decay-producing 
moment. For comparison the plus (+) symbols were computed from equation (7). The 
seven coefficients of the RFA of equation (7) were computed by substituting the 

A for the five values of k with C = 0 into equation (6). The RFA results com- 
22 

pare favorably with kernel- function results in this application. 

Figure 4(b) shows A^* weighted lift due to first-torsion-mode motion 

(note that the horizontal scale begins at 24). As with the moment A 2 2 # tll ^ s 

quantity displays more phase lag with decaying motion, and less phase lag and even 
a phase lead with growing motion in comparison to constant-amplitude motion. As in 
figure 4(a) the comparison between kernel-function and RFA results is good. 
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Flutter and damping-ratio boundaries .- Next the flutter boundary and selected 
damping-ratio boundaries for growing and decaying motion near the flutter boundary 
were calculated for this four-percent-thick rectangular model. This was done for 
four of the Mach numbers, namely, M = 0.756, 0.801, 0.856, and 0.904, listed in 
table III of Doggett, Rainey, and Morgan (1959). The model properties from their 
table 1 (a) for the first and second bending and first torsion modes were used; how- 
ever, because the modal damping coefficients g^ for first and second bending modes 
of the four-percent-thick wing are missing, those coefficients for the six-percent- 
thick wing were used from table 1(a). Calculated mode shapes for a uniform cantile- 
ver beam with midchord center of gravity and elastic axis were used. Moreover, the 
next two high frequency calculated modes, namely, second torsion and third bending, 
were included to ensure convergence of results with respect to the number of modes. 
For each of the four Mach numbers the density at the experimental flutter point was 
used. 


Figure 5 shows flutter boundary (*; = 0) as the flutter speed index C V/b o 03 a \|TI ) 
plotted against Mach number. The higher and lower boundaries are for £ = -0.025 
and 0.025, respectively. The comparison points from the RFA of equation (7) are very 
close to those from the kernel function. 


Two-Dimensional Wing Section 

Early studies of the DAST high-aspect-ratio transport wing model (see the fol- 
lowing section) produced unexpected wavy and looping curves of the complex general- 
ized forces as functions of reduced frequency for decaying motion. This behavior 
cast doubt on the adequacy of the computational process, including the convergence 
with respect to the downwash collocation order. Since collocation order can be 
extended very high for two-dimensional flow, aerodynamic forces were studied for the 
oscillating two-dimensional wing section to determine whether similar trends would be 
calculated and to test for convergence with respect to the number of downwash collo- 
cation points. Edwards (1979) presents similar examples of section lift due to 
plunging and pressure distributions resulting from both growing and decaying airfoil 
motions • 

For the present report the classical Possio integral equation formulation as 
described in Bland (1982) was used. The associated computer program for harmonic 
motion was generalized to growing and decaying exponential motion. For most of 
the two-dimensional section calculations, 64 downwash collocation points were 
used. Figure 6 shows the section lift due to plunging and pitching oscillations 
for M = 0.9 and for three motion-decay ratios, namely, (1) simple harmonic motion 
(£ = 0 ), (2) growing motion with 5 = - 0 . 1 , and (3) decaying motion with 5 = 0.1. 
Results from both the generalized kernel-function calculation (denoted ”K" ) and the 
RFA (based on 14 values of k) are included. The range of reduced frequency is 
large, from 0 to 4.0. 

Figure 6 (a) shows A^ , the generalized lift due to plunging. (The plunge 
half -amplitude is the semichord length.) The kernel-function results for £ = 0 
agree rather well with the RFA results of equation (7). The RFA results for grow- 
ing motion, £ = - 0 . 1 , also agree well with the kernel-function results. 

For decaying motion, £ = 0.1, the RFA of equation (7) is given by the smooth 
dashed curve. In great contrast is the looping curve of the kernel-function result. 
Certain values of k are labeled. As k increases beyond about 0.9, the curve 
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first develops a wavy behavior, then a near cusp, and finally widening loops. This 
unexpected wavy and looping behavior appears for values of k that are well beyond 
almost any wing flutter, but which might be reached for control-surface flutter. 

Figures 6(b) and 6(c) show the generalized lift A due to pitching about the 

quarter-chord with half -amplitude of 1 radian. Figure 6(b) shows the kernel-function 
result for £ = 0 and -0.1 from k = 0 to 5.0 and for £ = 0.1 from k = 0 to 0.5. 
It is apparent that the seven-term RFA function of equation (7) fitted to the 
14 values of this particular A^ ^ does not result in a good fit. This inadequate 

fit is due to the reversal of 5 ^ in the selected range of k and to the least- 

1 2 

square-error fit being over too great a range of k. An improvement would be needed 
for RFA applications. The spread for the three values of £ is as expected except 
for the looping near k = 0.5. Figure 6(c) gives an enlarged view of _ the behavior of 
A^ ^ f° r t ^ ie frequency range 0.5 < k < 4.5. The wavy behavior for £ = 0.1 begins 

near k = 0.5 and becomes looping behavior for k greater than about 1.0. 

As stated above, all the results shown were calculated with 64 collocation 
points on the chord. Checks were also made with 32 and 128 collocation points, and 
only insignificant differences were found even to k = 4.0. The conclusion is 
reached therefore that the kernel-function results shown are converged with respect 
to the number of collocation points, and that the looping behavior calculated for 
decaying motion with £ = 0.1 is a valid mathematical result. This looping behavior 
occurs for values of k well above the range of most wing flutter, except possibly 
control-surface flutter. It may be applicable to gust analysis that extends to high 
values of k. 


The DAST ARW-1 High-Aspect-Ratio Swept Wing 

The most realistic model analyzed here is the DAST ARW-1 (from Drones for Aero- 
dynamic and Structural Testing, aeroelastic research wing 1); see Murrow and Eckstrom 
(1979). The planform analyzed is shown in figure 1. It is essentially the planform 
shown in figure 2 of Newsom and Pototzky (1982), but no aerodynamic effects of hori- 
zontal or vertical tail surfaces or of a fuselage are included. The 12 spanwise 
symmetric mode shapes and frequencies employed are those calculated for the original 
DAST ARW-1 wing via the NASTRAN 2 finite-element analysis that was reported in Newsom 
and Pototzky (1982). These modes consist of 2 rigid-body (vertical translation and 
pitch) and the first 10 elastic modes, not including any aileron motion. The aileron 
hinge is treated as locked (no aileron rotation). 

Aerodynamic-force results .- Figure 7 shows the A^ from kernel-function 
results and the RFA (based on nine values of k) for the damping ratios £ = 0, ±0.1, 
and ±0.577 (tan(±30°)) for 0 < k < 1.0 and M = 0.8. Mode 5 (the third elastic 
mode), although strongly coupled, is judged to be best described as "first torsion." 
The original ARW-1 wing described by Edwards (1979) and Newsom and Pototzky (1982) 
experienced flutter during flight testing at k * 0.16. 


^NASTRAN is a registered trademark of the National Aeronautics and Space 
Administration • 
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The looping behavior of the kernel-function result for 5 « 0.577 is like the 
looping of the two-dimensional case described in a preceding section and is believed 
to be a valid result. The RFA result agrees with the kernel-function result except 
for the highly damped motion (£ = 0.577) at values of k about 0.2 and higher, 
especially for k higher than 0.3. 

Flutter analysis results .- The present flutter analysis employs the subsonic 
kernel-function aerodynamics generalized to growing and decaying motion. Newsom 
and Pototzky (1982) reported the flutter and preflutter results obtained based on 
doublet-lattice aerodynamics for simple harmonic motion and extrapolation into the 
complex s-plane via the RFA. 

The technique of strengthening the analytical aerodynamic forces by multiplying 

by the ratio of wind-tunnel static lift to the analytical static lift employed by 

Newsom and Pototzky (1982) is used in the present analysis. The upper curve of 

C versus M in figure 8 herein is from figure 17 of Byrdsong and Hallissy (1979) 

1 * 

a 

and is for the complete model including the horizontal tail. Since the analysis here 
is of the wing only, the second curve of figure 8 was obtained from the tail-off data 
(plots of C L vs. a) of figures 11(a) to 11(e) of Byrdsong and Hallissy (1979). 

The lowest curve of figure 8 is from the kernel-function program using program- 
default downwash collocation at 6 chord stations on each of 12 span stations. The 
difference in percent from the bottom (analytical) curve up to the tail-off experi- 
mental curve is indicated for several Mach numbers; this difference was used to 
strengthen the generalized force matrix in the flutter analysis. 

Flutter stability was analyzed at Mach numbers from 0.70 up to 0.82 for 
15000 feet and up to 0.94 for 25000 feet. Densities and sound speeds were used from 
the U.S. Standard Atmosphere, 1962. Figure 9(a) for 15000 feet and figure 9(b) for 
25000 feet show the calculated damping ratio 5 and associated frequency f ver- 
sus Mach number. The circles are the experimental values that appear in figures 8 
and 10 of Newsom and Pototzky (1982), and the flight flutter point from their fig- 
ure 10 is shown by the arrow. The good correlation of the present analysis and 
experiment is evident. The analytical results of Newsom and Pototzky (1982), based 
on doublet-lattice aerodynamics with the strengthening of their figure 5, are given 
by the dashed curves. As can be observed, these results are comparable to those of 
the present analysis and to the flight experiments. 


CONCLUDING REMARKS 

The subsonic kernel function has been generalized to the s-plane for growing and 
decaying oscillatory motion. The function was substituted into a flutter analysis 
program, and that program was also adapted throughout to obtain the desired solutions 
for constant-damping-ratio boundaries, including the usual zero-decay flutter bound- 
ary. A rational function approximation (RFA) used in control theory for approximat- 
ing s-plane aerodynamic forces was used for some comparisons. 

For the cases analyzed, the RFA aerodynamics gave (or could give) very good 
agreement with the s-plane subsonic-kernel-function results for the lower reduced 
frequencies that characterize wing flutter. But for combinations of higher damping 
ratios and of higher frequencies that often characterize control-surface flutter, the 
RFA aerodynamics did not follow the looping and spiraling kernel -function results. 
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For a clipped-tip delta -wing flutter model, the s-plane kernel function pre- 
dicted a less steep gradient of damping ratio versus density as the flutter boundary 
was approached and exceeded than did the conventional V-g solution method that uses 
simple harmonic aerodynamics. 

An aspect-ratio-five rectangular-planform model was analyzed for four Mach num- 
bers (M) ranging from 0.756 to 0.904. Constant-damping-ratio boundaries, one for 
decaying and one for growing motion, are presented near the flutter boundary. At 
M = 0.904, study of the weighted torsional moment due to first-tors ion-mode oscilla- 
tion showed that in comparison to the phase lagging moment of constant-amplitude 
motion, decay produced a greater lag of that moment, while growth produced a lesser 
lag. Thus, in terms of this first torsion moment only, decay promotes faster decay, 
and growth promotes slower decay; but for coupled-mode flutter there is no assurance 
that the overall flutter motion will show this same effect. The RFA aerodynamics 
gave good comparisons for both aerodynamic forces and decay-ratio results for the low 
reduced frequencies studied. 

3 

An analysis of the DAST ARW-1 transport-type wxng was also made. Plots of 
damping ratio versus flight Mach number were calculated for two altitudes, 15000 and 
25000 feet. Good comparisons of damping ratios and associated frequencies were 
obtained with flight experiments and also with a doublet-lattice analysis that used 
the RFA. The same type of analytical treatment was used here as was used in a refer- 
enced doublet-lattice analysis; namely, in the flutter analysis the generalized 
forces were strengthened by multiplying by the ratio of the lift-curve slope C 

Xj 

a 

obtained on a wind-tunnel model (tail off) to that from the subsonic-kernel-function 
aerodynamics. For M = 0.8 and reduced frequencies ranging from 0 to 1.0 the gener- 
alized w f irs t- tors ion-mode M forces calculated from the kernel function agreed with 
those calculated from the rational function approximation (as in Abel, NASA TP-1367, 
1979, for example) except at large damping ratios. 

At reduced frequencies higher than that of the experimental flutter, an unex- 
pected doubling back and looping of some of the generalized forces for strongly 
decaying motion were found in the complex plane. This behavior cast doubt on the 
adequacy of the computational process for the DAST wing, including the convergence 
with respect to the downwash collocation order. Therefore the two-dimensional-flow 
wing section was studied because the collocation order can be extended very high 
and collocation convergence can be studied. The classical Possio integral equation 
formulation was used. An existing computer program for harmonic motion was general- 
ized to growing and decaying motion. For moderately decaying motion a looping behav- 
ior of the lift forces due to both translational and pitching oscillations was found 
at higher reduced frequencies for collocation-converged results. Because of these 
two-dimensional-flow results the similar behavior for the DAST ARW-1 in three- 
dimensional flow at high reduced frequencies and high damping ratio is concluded to 
be valid. In any event the reduced frequency of flutter is well below that of the 
looping behavior. The RFA aerodynamics did not follow the looping behavior at higher 
reduced frequencies. 

In the appendix, procedures were described and evaluated for the accurate calcu- 
lation of the kernel function for general s-plane motion. A 12-term exponential 


3 

wing 1 


Drones for Aerodynamic and Structural Testing (DAST) , 
(ARW-1 ). 


aeroelastic research 


13 



series approximation in the kernel-function integral was shown to provide good accu- 
racy for growing and decaying motion in the s-plane extending at least 45° on both 
sides of the imaginary s-axis that represents harmonic motion. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton , VA 23665 
February 17, 1984 
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APPENDIX 


CALCULATION OF THE SUBSONIC KERNEL FUNCTION IN THE s-PLANE 

This appendix serves four purposes* It describes the algorithm used in this 
report to compute the kernel of the downwash integral equation. It indicates in what 
part of the complex s-plane the algorithm gives acceptable results. It explains why 
this particular algorithm was used. It also explains why the algorithm used, or any 
similar algorithm, cannot be expected to give good results for all complex values of 
reduced frequency. 

The kernel function, essentially as in equation (3) but with k replaced by 
k c , is 


K(M,k ,X ,Y > 

K(M,k ,x ,y ) = r — - 

coo n 2 2 

* 


where 


K(M,k ,x ,y ) = exp(-ik x ) 
coo co 


1 + ^Jexp(-ik c u 1 |y o |) 


I, 


- ik c' y o'J ( 1 " J" j 

'u \ M 1 + u 


u 


exp(-ik |y |u) du 


c ' o' 


(A1 ) 


and where k c is the complex reduced frequency 


k c “ k(1 + i? ) 


and 


. 2 2 2 
R = \ x + 0 y 
\ o o 


u i = (_x o + MR) A B2 ly 0 ^ 


In equation (A1 ) let 


s = ik^ = si /V 


(A2 ) 
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Then 


K(M,-is,x ,y ) = exp ( -sx ) 
o o o 


1 + exp(-su 1 |y o |) - p(5|y o l) 


(A3) 


where 


F(z) 




1 - 


\ 1 + 


exp(-zu) du 


(A4) 


and the symbol z = s|y Q | is not the z of the main text. Integral equation (A4) 
defines the function F(z) only for | arg z| < ir/2. For the multi leaved Riemann 
surface, (arg z| > tt/2, F(z) is defined by analytic continuation as follows: 


where 


F(z) = F (z) - F (z) 


r (z) =z I (1 - - r ) exp(-zu) 
1 J 0 \ >|l + u / 


du 


u 


F^ ( z) = z 


f ' 

U _ . 

- u 1 

K 

[ ' 

0 + u 2 / 


exp(-zu) du 


By integrating by parts and applying equations 12.1.8, 12.1.9, and 9.1.5 of 
Abramowitz and Stegun (1964), one obtains 


F^(z) = 1 + z - ~ z ( z ) + z Y^(z) 


(A5) 


where H^(z) is a Struve function and Y-j (z) is a Bessel function of the second 
kind. To show the analytic behavior of F 2 (z), expand exp(-zu) into a power series 
and integrate termwise to obtain 


-E 


, .m 

(-1) g 


F 2 (Z) ^ .! 

m=0 


m m+1 
z 


(A6) 
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where 



Since |g m | is bounded by |u*|| m+1 /(m + 1), the series (eq. (A6)) has an infinite 
radius of convergence , and hence F 2 (z) an ent ^ re analytic function# In equa- 

tion (A5), the Struve function z H*j(z) is also an entire analytic function. The 
Bessel function z (z) can be expressed 

z Y (z) = — z J (z) In + G(z) 

1 tt 1 2 

where z J-| (z) and G(z) are both entire analytic functions of z. Thus equa- 
tions (A5 ) and (A6) furnish an analytic continuation of F(z) into the multi leaved 
Riemann surface | arg z| > tt/2. 

There are no singularities in F(z) in the finite part of the complex plane 
except for a logarithmic branch point at z = 0. The natural branch cut of F(z) 
is the natural branch cut of In z/2, namely , the negative real z-axis. 

The integral of equation (A4) is evaluated numerically by replacing the alge- 
braic part of the integrand by an exponential approximation and integrating termwise: 


1 


u 


1 + u 



exp (-2 


V) 


(0 < u < °°) 


(A7) 


An approximation of this sort is used for several reasons: (1) it is a function of 

a single variable; (2) only one exponential is needed per kernel function, contribut- 
ing to economy of calculation? (3) a modest number of terms provide nearly single- 
precision computer-word accuracy for harmonic motion? and (4) as described below, 
good accuracy extends to rather large damping ratios. The coefficients in equa- 
tion (A7), namely. 


b = 0.009054814793 


a 1 

= 

0.000319759140 

a 2 

= 

-0.000055461471 

4 

= 

0.002726074362 

a 4 


0.005749551566 

a 5 

= 

0.031455895072 

4 

= 

0.106031126212 

a 7 

= 

0.40683801 1567 

a 8 

- 

0.798112357155 

a 9 

= 

-0.417749229098 

a 10 

= 

0.077480713894 

a 11 


-0.012677284771 

a 1 2 

= 

0.001787032960 
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are for the approximation (eq. (A7)) designated D12.1 by Desmarais (1982), Then, if 

u 1 > 0, 


p(i|y o l) “ s | y Q | 






2 b + s I y 


(A8) 


where 


e 1 = exp(-bu 1 ) 

e i E ( ®A-1 )2 (1=2 to 12) 


The closed-form integrals that were used to derive equation (A8) all converge only if 
Re(s|y 0 |) > -b • For lower values of Re(s|y 0 |)/ equation (A8) is deduced by analytic 
continuation. 

Since approximation D12.1 is valid only for u > 0, another approximation to 
F(s|y 0 |) is used for u^ < 0: 


12 / 12 
p(i|y |) - -2 - 2, ;y /V * 4- «p(-So |y I) (2 + ^ 

(2 b) - (sy^) \ f 




1-1 2 b - s|y o' 


(A9) 


where 

e 1 = exp(bu 1 ) 

= (e A _ J| ) 2 (1 = 2 to 12) 

Approximation (A8) has poles at 

sly I = -2^b ( £ = 1 to 12) 

'o' 

that is, on the negative real s-axis, and approximation (A9) has poles at 

sly | = ±2 £ b (l = 1 to 12) 

1 o 1 
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that is, on the negative and positive real s-axes. The function F(s|y 0 |) has no 
singularities except a branch point at s|y Q | = 0. Thus one would expect approx- 
imation (A8 ) to become unusable as arg s approaches ±ir, and approximation (A9) 
to behave similarly as arg 5 approaches either ±tt or 0. This is shown in fig- 
ure 10(a). These are plots of the real and imaginary parts of K versus |k c y 0 | 
for |k c y 0 | = 0 to 8 and u-j = 0. The exact ic (solid line) is compared with the 
approximation (dashed line). The value of u^ is zero so that the exact K is 
easily computed from equation (A5). For 0 = tt/ 8, the exact and approximate curves 
are almost indistinguishable. For 0 = ir/4, the error is noticeable but approxima- 
tion (A7) gives acceptable engineering accuracy. Figure 10(b) shows the same infor- 
mation for an earlier and commonly used exponential approximation similar to D12.1 
presented in Laschka (1970) and designated L1 1 in Desmarais (1982). Approxima- 
tion L1 1 is about as bad at 0 = tt/8 as approximation D12.1 is at 0 = tt/4. At 
0 = tt/4, approximation L1 1 is unacceptable. The lower plot of figure 10(b) also 
appears as part of figure 1 of Ashley and Boyd (1980) along with a still earlier 
approximation from Watkins, Woolston, and Cunningham (1959) designated W4 in 
Desmarais (1982). In Ashley and Boyd (1980), the poor performance of L1 1 and W4 for 
0 = tt/4 is attributed to the proximity of the poles of the computed kernel. How- 
ever, approximation D12.1, which has poles along the same line as L1 1 , performs very 
well at 0 = tt/4. This is so because the error induced near a pole is proportional 
to the magnitude of the residue at that pole, and each residue is proportional to 
the associated coefficient a^ . For approximation L1 1 , maxla^l = 644.8, whereas for 
approximation D12.1, maxla^) = 0.798. 

To sum up, if 1 0 1 < tt/4, then equations (A8) and (A9), which use approximation 
D12.1, give acceptable accuracy. All calculations in this report used values of k c 
for which 1 0 [ was much less than tt/4 and were performed using approximation D12.1 
except for a few check calculations that were performed using approximation D24.2 of 
Desmarais (1982); and no significant differences were found. Approximation D24.2 is 
more accurate than D12.1 and takes about twice as long to execute. 

It is the authors' opinion that exponential approximations, such as D12.1, for 
which max|a k | < 1, provide a satisfactory way of computing K(M,k c ,x Q ,y Q ) (where 

k c = -is), to the accuracy needed for aerodynamic-f orce calculations, in the half 
of the complex s-plane for which either | —0 1 < tt/4 or ( 9 1 < tt/4; that is, for 
tt / 4 < arg s < 3tt/4 and -3tt/4 < arg s < -ir/4. For the remainder of the complex 
s-plane, | arg s| < tt/4 or |arg s| > 3tt/4, estimates of K based on exponential 

approximations to 1 - (u/\)1 + u 2 ) become inaccurate or uneconomical, or both. 
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Figure 4.- Generalized aerodynamic forces for growing and decaying motions of the 
four-percent-thick rectangular-planform flutter model of Doggett, Rainey, and 
Morgan (NASA TM X-79, 1959). 0 < k < 0.228; M = 0.904. 










(a) A , lift due to plunging oscillation. 0 < k < 4. 
11 

Figure 6.- Generalized aerodynamic forces for growing and 
decaying motions of a two-dimensional airfoil section. 






(c) A for 0*5 < k < 4,5 with additional kernel-function results 
12 for Z = 0.1. 

Figure 6*- Concluded* 
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for growing and decaying motions of the 


DAST ARW-1 wing. M = 0.8; 0 < k < 1.0. Mode 5 is "first torsion 
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(b) Altitude = 25000 feet. 
Figure 9.- Concluded. 
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